#####################################
### Matthew Cebul + Sharan Grewal ###
### Conscription & NVCs           ###
### Analysis                      ###
#####################################

rm(list = ls())
library(plyr)
library(tidyverse)
library(reshape2)
library(stargazer)
library(ggplot2)
library(sandwich)
library(margins)
library(prediction)
require(survival)

setwd("/Users/mcebul/Dropbox/Conscription/Submissions/CPS/Final")

# Load Data
data <- read_rds("LargeN_clean.rds")

# Drop country-years where regime lacks a standing army
dat <- data[which(data$mil==1),] 

# set WD for storing images 
setwd("/Users/mcebul/Dropbox/Conscription/Submissions/CPS/Final/Figures")


##############
## Analysis ##
##############

######
# TABLE 1, NVONSET
######

# Nonviolent Campaign Onset, Models
NVOnset1 <- lm(NVonset~conscript, data=dat)
NVOnset2 <- lm(NVonset~conscript+ # conscripted army (Toronto)
                 War+             # Interstate War Indicator (Reiter, Stam, Horowitz 2016)
                 Rivalry+         # Interstate Rivalry Indicator (Thompson Dreyer 2012)
                 milperper+       # mil personnel (COW_NMC)
                 logmilex2+       # mil spending (COW_NMC)
                 demo_r+          # Democracy dummy (Magaloni)
                 milregime+       # Military regime type dummy (Magaloni)
                 BritCol+         # British Colonial History (Asal + Toronto)
                 gdploglag+       # lagged + logged GDPPC (PWT + Maddison)
                 chgdpcomb+       # Change GDPPC (PWT + Maddison)
                 tpoplog+         # Population (COW_NMC)
                 youth+           # 15-24 as % of total (WDI)
                 urban+           # % Urban (COW_NMC)
                 Mobilelog+       # Mobile phone per 100 (WDI)
                 diffuselog+      # Nonviolent onsets in region (NAVCO)
                 Past_success+    # Binary indicator of past successful NV campaign (NAVCO)
                 Past_defect+      # Binary indicator of military defections against past campaigns (NAVCO)
                 factor(region)+  # Regional indicator
                 factor(year),    # Fixed effect for year
               data=dat)

Models_ses <- lapply(list(NVOnset1, NVOnset2), function (x) sqrt(diag(vcovHC(x, type = "HC1"))))
labels <- c("Conscript", "War", "Rivalry", "MilSize (log)", "MilSpend (log)", "Democracy", "MilRegime", 
            "BritCol", "GDP (log)", "GDPChange", "PopSize", "Youth", "Urban", "Mobile (log)", "Diffusion (log)",
            "PastSuccess", "PastDefection", "Region(Americas)", "Region(East Asia \\& Pacific)", "Region(Europe \\& Eurasia)", 
            "Region(Middle East \\& North Africa)", "Region(South \\& Central Asia)")

stargazer(NVOnset1, NVOnset2, title ="Conscription and Nonviolent Campaigns", 
          se = Models_ses,
          dep.var.caption = "Dependent Variable",
          dep.var.labels = c("Nonviolent Campaign Onset"),
          align=TRUE, float = TRUE, column.sep.width = "0pt", single.row = TRUE,
          #order=paste0("^", vars.order , "$"),
          covariate.labels = labels,
          omit.stat = c("f", "ser"),
          notes = c("Robust SEs computed via Huber-White sandwich estimator."))


#####
# FIGURE 2 (PREDICTED VALUES NVONSET)
##### 
pred1 <- prediction(NVOnset2, at= list(conscript = c("0","1")), level = 0.95)
print(pred1sum <- summary(pred1))
colnames(pred1sum)[1] <- "conscript"

#Plot it
pdf("NVOnset_Predicted.pdf")
ggplot(pred1sum, aes(x = conscript, y = Prediction)) +
geom_col(aes(x = conscript), fill = "royalblue2", alpha = 0.75) +
  geom_errorbar(aes(x = conscript, ymin = lower,
                     ymax = upper, width = 0.1)) +
  ylim(0,0.035) +
  xlab("Conscription") + ylab("Prob. NVOnset") +
  geom_text(aes(x=conscript, y = Prediction,
                label=round(Prediction, 4)), size=5, hjust = 1.25, vjust=-.5) +
  theme_bw(base_size=16)
dev.off()


#####
# TABLE 2, SIZE DEFECT SUCCESS
#####

# Nonviolent Campaign Size
NVSize2 <- lm(Size2~conscript+ # conscripted army (Toronto)
                War+             # Interstate War Indicator (Reiter, Stam, Horowitz 2016)
                Rivalry+         # Interstate Rivalry Indicator (Thompson Dreyer 2012)
                milperper+       # mil personnel (COW_NMC)
                logmilex2+       # mil spending (COW_NMC)
                demo_r+          # Democracy dummy (Magaloni)
                milregime+       # Military regime type dummy (Magaloni)
                BritCol+         # British Colonial History (Asal + Toronto)
                gdploglag+       # lagged + logged GDPPC (PWT + Maddison)
                chgdpcomb+       # Change GDPPC (PWT + Maddison)
                tpoplog+         # Population (COW_NMC)
                youth+           # 15-24 as % of total (WDI)
                urban+           # % Urban (COW_NMC)
                Mobilelog+       # Mobile phone per 100 (WDI)
                diffuselog+      # Nonviolent onsets in region (NAVCO)
                Past_success+    # Binary indicator of past successful NV campaign (NAVCO)
                Past_defect+     # Binary indicator of military defections against past campaigns (NAVCO)
                factor(region)+  # Regional indicator
                factor(year),    # Fixed effect for year
              data=dat[dat$prim_meth==1,])
Defect2 <- lm(defect~conscript+
                War+             # Interstate War Indicator (Reiter, Stam, Horowitz 2016)
                Rivalry+         # Interstate Rivalry Indicator (Thompson Dreyer 2012)
                milperper+       # mil personnel (COW_NMC)
                logmilex2+       # mil spending (COW_NMC)
                demo_r+          # Democracy dummy (Magaloni)
                milregime+       # Military regime type dummy (Magaloni)
                BritCol+         # British Colonial History (Asal + Toronto)
                gdploglag+       # lagged + logged GDPPC (PWT + Maddison)
                chgdpcomb+       # Change GDPPC (PWT + Maddison)
                tpoplog+         # Population (COW_NMC)
                youth+           # 15-24 as % of total (WDI)
                urban+           # % Urban (COW_NMC)
                Mobilelog+       # Mobile phone per 100 (WDI)
                diffuselog+      # Nonviolent onsets in region (NAVCO)
                Past_success+    # Binary indicator of past successful NV campaign (NAVCO)
                Past_defect+     # Binary indicator of military defections against past campaigns (NAVCO)
                factor(region)+  # Regional indicator
                factor(year),    # Fixed effect for year
              data=dat[(dat$camp_goals>= 0 & dat$camp_goals <=4) & dat$prim_meth==1,])
Success2 <- lm(success~conscript+
                 War+             # Interstate War Indicator (Reiter, Stam, Horowitz 2016)
                 Rivalry+         # Interstate Rivalry Indicator (Thompson Dreyer 2012)
                 milperper+       # mil personnel (COW_NMC)
                 logmilex2+       # mil spending (COW_NMC)
                 demo_r+          # Democracy dummy (Magaloni)
                 milregime+       # Military regime type dummy (Magaloni)
                 BritCol+         # British Colonial History (Asal + Toronto)
                 gdploglag+       # lagged + logged GDPPC (PWT + Maddison)
                 chgdpcomb+       # Change GDPPC (PWT + Maddison)
                 tpoplog+         # Population (COW_NMC)
                 youth+           # 15-24 as % of total (WDI)
                 urban+           # % Urban (COW_NMC)
                 Mobilelog+       # Mobile phone per 100 (WDI)
                 diffuselog+      # Nonviolent onsets in region (NAVCO)
                 Past_success+    # Binary indicator of past successful NV campaign (NAVCO)
                 Past_defect+     # Binary indicator of military defections against past campaigns (NAVCO)
                 factor(region)+  # Regional indicator
                 factor(year),    # Fixed effect for year
               data=dat[(dat$camp_goals>= 0 & dat$camp_goals <=4) & dat$prim_meth==1,])

# Make Table
Models_ses <- lapply(list(NVSize2, Defect2, Success2), function (x) sqrt(diag(vcovHC(x, type = "HC1"))))
stargazer(NVSize2, Defect2, Success2, title ="Campaign-Year Analysis", 
          se = Models_ses,
          dep.var.caption = "Dependent Variable",
          dep.var.labels = c("NV Campaign Size", "Defection", "NV Campaign Success"),
          covariate.labels = c("Conscript", "Constant"),
          align=TRUE, float = TRUE, column.sep.width = "0pt",
          keep=c("conscript1", "Constant"),
          omit.stat = c("f", "ser"),
          notes = c("Robust SEs computed via Huber-White sandwich estimator."))


######
# TABLE 3, VIOLENCE ONSET + TACTICS
######

VioOnset2 <- lm(VioOnset~conscript+ # conscripted army (Toronto)
                  War+             # Interstate War Indicator (Reiter, Stam, Horowitz 2016)
                  Rivalry+         # Interstate Rivalry Indicator (Thompson Dreyer 2012)
                  milperper+       # mil personnel (COW_NMC)
                  logmilex+        # mil spending (COW_NMC)
                  demo_r+          # Democracy dummy (Magaloni)
                  milregime+       # Military regime type dummy (Magaloni)
                  BritCol+         # British Colonial History (Asal + Toronto)
                  gdploglag+       # lagged + logged GDPPC (PWT + Maddison)
                  chgdpcomb+       # Change GDPPC (PWT + Maddison)
                  tpoplog+         # Population (COW_NMC)
                  youth+           # 15-24 as % of total (WDI)
                  urban+           # % Urban (COW_NMC)
                  Mobilelog+       # Mobile phone per 100 (WDI)
                  diffuselog+      # Nonviolent onsets in region (NAVCO)
                  Past_success+    # Binary indicator of military defections against past campaigns (NAVCO)
                  Past_defect+     # Binary indicator of military defections against past campaigns (NAVCO)
                  factor(region)+  # Regional indicator
                  factor(year),    # Year Fixed effect,  
                data=dat)
Tactics2 <- lm(prim_meth~conscript+
                 War+             # Interstate War Indicator (Reiter, Stam, Horowitz 2016)
                 Rivalry+         # Interstate Rivalry Indicator (Thompson Dreyer 2012)
                 milperper+       # mil personnel (COW_NMC)
                 logmilex+        # mil spending (COW_NMC)
                 demo_r+          # Democracy dummy (Magaloni)
                 milregime+       # Military regime type dummy (Magaloni)
                 BritCol+         # British Colonial History (Asal + Toronto)
                 gdploglag+       # lagged + logged GDPPC (PWT + Maddison)
                 chgdpcomb+       # Change GDPPC (PWT + Maddison)
                 tpoplog+         # Population (COW_NMC)
                 youth+           # 15-24 as % of total (WDI)
                 urban+           # % Urban (COW_NMC)
                 Mobilelog+       # Mobile phone per 100 (WDI)
                 diffuselog+      # Nonviolent onsets in region (NAVCO)
                 Past_success+    # Binary indicator of military defections against past campaigns (NAVCO)
                 Past_defect+  # Binary indicator of military defections against past campaigns (NAVCO)
                 factor(region)+  # Regional indicator
                 factor(year),    # Year Fixed effect,  
               data=dat[(dat$camp_goals>= 0 & dat$camp_goals <=4),])

#Make table
Models_ses <- lapply(list(VioOnset2, Tactics2), function (x) sqrt(diag(vcovHC(x, type = "HC1"))))
stargazer(VioOnset2, Tactics2, title ="Nonviolence vs. Violence", 
          se = Models_ses,
          dep.var.caption = "Dependent Variable",
          dep.var.labels = c("Vio Campaign Onset", "Campaign Tactics"),
          covariate.labels = c("Conscript", "Constant"),
          align=TRUE, float = TRUE, column.sep.width = "0pt",
          keep=c("conscript1", "Constant"),
          omit.stat = c("f", "ser"),
          notes = c("Robust SEs computed via Huber-White sandwich estimator."))



#############
## APPENDIX
#############


#####
# Full Table, Campaign-Level Analysis
#####

Models_ses <- lapply(list(NVSize2, Defect2, Success2), function (x) sqrt(diag(vcovHC(x, type = "HC1"))))
labels <- c("Conscript", "War", "Rivalry", "MilSize (log)", "MilSpend (log)", "Democracy", "MilRegime", 
            "BritCol", "GDP (log)", "GDPChange", "PopSize", "Youth", "Urban", "Mobile (log)", "Diffusion (log)",
            "PastSuccess", "PastDefection", "Region(Americas)", "Region(East Asia \\& Pacific)", "Region(Europe \\& Eurasia)", 
            "Region(Middle East \\& North Africa)", "Region(South \\& Central Asia)")

stargazer(NVSize2, Defect2, Success2, title ="Campaign-Year Analysis, Full", 
          se = Models_ses,
          dep.var.caption = "Dependent Variable",
          dep.var.labels = c("NV Campaign Size", "Defection", "NV Campaign Success"),
          align=TRUE, float = TRUE, column.sep.width = "0pt", single.row = TRUE,
          covariate.labels = labels,
          omit.stat = c("f", "ser"),
          notes = c("Robust SEs computed via Huber-White sandwich estimator."))


#####
# Full Table, Vio + Tactics Analysis
#####

Models_ses <- lapply(list(VioOnset2, Tactics2), function (x) sqrt(diag(vcovHC(x, type = "HC1"))))
labels <- c("Conscript", "War", "Rivalry", "MilSize (log)", "MilSpend (log)", "Democracy", "MilRegime", 
            "BritCol", "GDP (log)", "GDPChange", "PopSize", "Youth", "Urban", "Mobile (log)", "Diffusion (log)",
            "PastSuccess", "PastDefection", "Region(Americas)", "Region(East Asia \\& Pacific)", "Region(Europe \\& Eurasia)", 
            "Region(Middle East \\& North Africa)", "Region(South \\& Central Asia)")

stargazer(VioOnset2, Tactics2, title ="NV vs. Vio Tactics", 
          se = Models_ses,
          dep.var.caption = "Dependent Variable",
          dep.var.labels = c("Vio Campaign Onset", "Campaign Tactics"),
          covariate.labels = labels,
          align=TRUE, float = TRUE, column.sep.width = "0pt", single.row = TRUE,
          omit.stat = c("f", "ser"),
          notes = c("Robust SEs computed via Huber-White sandwich estimator."))


#####
# NVONSET EXTENSION + NV COUNTS
#####

# NVONSET
NVOnset3 <- lm(NVonset~conscript+ # conscripted army (Toronto)
                 War+             # Interstate War Indicator (Reiter, Stam, Horowitz 2016)
                 Rivalry+         # Interstate Rivalry Indicator (Thompson Dreyer 2012)
                 milperper+       # mil personnel (COW_NMC)
                 logmilex2+        # mil spending (COW_NMC)
                 demo_r+          # Democracy dummy (Magaloni)
                 milregime+       # Military regime type dummy (Magaloni)
                 BritCol+         # British Colonial History (Asal + Toronto)
                 gdploglag+       # lagged + logged GDPPC (PWT + Maddison)
                 chgdpcomb+       # Change GDPPC (PWT + Maddison)
                 tpoplog+         # Population (COW_NMC)
                 youth+           # 15-24 as % of total (WDI)
                 urban+           # % Urban (COW_NMC)
                 Mobilelog+       # Mobile phone per 100 (WDI)
                 diffuselog+      # Nonviolent onsets in region (NAVCO)
                 Past_success+    # Binary indicator of past successful NV campaign (NAVCO)
                 Past_defect+     # Binary indicator of military defections against past campaigns (NAVCO)
                 spons+           # Indicator for authoritarian client regimes (CASEY)
                 PHYSINT+         # CIRI Physical Integrity Index
                 factor(region)+  # Regional indicator
                 factor(year),    # Fixed effect for year
               data=dat)

# COUNTS
dat$countryyear <- paste(dat$country, dat$year) #prepare to drop duplicate country-years
NVCounts <- lm(NVOnsetSum~conscript+ # conscripted army (Toronto)
                 War+             # Interstate War Indicator (Reiter, Stam, Horowitz 2016)
                 Rivalry+         # Interstate Rivalry Indicator (Thompson Dreyer 2012)
                 milperper+       # mil personnel (COW_NMC)
                 logmilex2+       # mil spending (COW_NMC)
                 demo_r+          # Democracy dummy (Magaloni)
                 milregime+       # Military regime type dummy (Magaloni)
                 BritCol+         # British Colonial History (Asal + Toronto)
                 gdploglag+       # lagged + logged GDPPC (PWT + Maddison)
                 chgdpcomb+       # Change GDPPC (PWT + Maddison)
                 tpoplog+         # Population (COW_NMC)
                 youth+           # 15-24 as % of total (WDI)
                 urban+           # % Urban (COW_NMC)
                 Mobilelog+       # Mobile phone per 100 (WDI)
                 diffuselog+      # Nonviolent onsets in region (NAVCO)
                 Past_success+    # Binary indicator of past successful NV campaign (NAVCO)
                 Past_defect+     # Binary indicator of military defections against past campaigns (NAVCO)
                 spons+           # Indicator for authoritarian client regimes (CASEY)
                 PHYSINT+         # CIRI Physical Integrity Index
                 factor(region)+  # Regional indicator
                 factor(year),    # Fixed effect for year
               data=dat[!duplicated(dat$countryyear),])

Models_ses <- lapply(list(NVOnset3, NVCounts), function (x) sqrt(diag(vcovHC(x, type = "HC2"))))
labels <- c("Conscript", "War", "Rivalry", "MilSize (log)", "MilSpend (log)", "Democracy", "MilRegime", 
            "BritCol", "GDP (log)", "GDPChange", "PopSize", "Youth", "Urban", "Mobile (log)", "Diffusion (log)",
            "PastSuccess", "PastDefection", "ForeignSponsor", "PhysInt", "Region(Americas)", "Region(East Asia \\& Pacific)", "Region(Europe \\& Eurasia)", 
            "Region(Middle East \\& North Africa)", "Region(South \\& Central Asia)")

stargazer(NVOnset3, NVCounts, title ="NVOnset Robustness Checks", 
          se = Models_ses,
          dep.var.caption = "Dependent Variable",
          dep.var.labels = c("NV Campaign Onset", "NVOnset Count"),
          align=TRUE, float = TRUE, column.sep.width = "0pt", single.row = TRUE,
          covariate.labels = labels,
          omit.stat = c("f", "ser"),
          notes = c("Robust SEs computed via Huber-White sandwich estimator."))


#####
# Duration Analysis, Campaign Level
#####
 
# First, a coarse categorization of campaigns as either nonviolent or violent
# if more campaign years nonviolent than violent, we code as nonviolent.
dat <- dat %>%
  group_by(camp_name) %>%
  mutate(type = ifelse(mean(prim_meth) >= 0.5, 1,0)) #type=1=predominantly nonviolent

# Now an indicator for whether the campaign end-date is right censored.
# NAVCO 2.1 provides this in "status", which indicates if campaign was ongoing
# after the end of their observation period (0 = campaign ended, 1 = campaign ongoing). 
# The Surv() model object takes status as 0 = alive/no event, 1 = dead/event. 
# So we need to flip status.
dat$status2 <- ifelse(dat$status==1, 0, 1)

# Now, trim dataset to one row per NAVCO campaign, keeping only the first campaign-year
dataslim <- dat %>% 
  group_by(camp_name) %>%
  filter(row_number()==1)

# Cox Proportional Hazards Models, Nonviolent Campaigns
Duration1 <- coxph(Surv(Duration, status2) ~ conscript, data=dataslim[dataslim$type==1,])
Duration2 <- coxph(Surv(Duration, status2) ~ conscript+ 
                     War+             # Interstate War Indicator (Reiter, Stam, Horowitz 2016)
                     Rivalry+         # Interstate Rivalry Indicator (Thompson Dreyer 2012)
                     milperper+       # mil personnel (COW_NMC)
                     logmilex+        # mil spending (COW_NMC)
                     demo_r+          # Democracy dummy (Magaloni)
                     milregime+       # Military regime type dummy (Magaloni)
                     BritCol+         # British Colonial History (Asal + Toronto)
                     gdploglag+       # lagged + logged GDPPC (PWT + Maddison)
                     chgdpcomb+       # Change GDPPC (PWT + Maddison)
                     tpoplog+         # Population (COW_NMC)
                     youth+           # 15-24 as % of total (WDI)
                     urban+            # % Urban (COW_NMC)
                     Mobilelog+       # Mobile phone per 100 (WDI)
                     diffuselog+      # Nonviolent onsets in region (NAVCO)
                     Past_success+    # Binary indicator of military defections against past campaigns (NAVCO)
                     Past_defect+     # Binary indicator of military defections against past campaigns (NAVCO)
                     factor(region),  # Regional indicator
                     #factor(year),    # Year Fixed effect, dropped b/c limited degrees of freedom
                     data=dataslim[dataslim$type==1,])


DurationTable <- summary(Duration2)$coefficients
rownames(DurationTable) <- c("Conscript", "War", "Rivalry", "MilSize (log)", "MilSpend (log)", "Democracy", "MilRegime", 
                             "BritCol", "GDP (log)", "GDPChange", "PopSize", "Youth", "Urban", "Mobile (log)", "Diffusion (log)",
                             "PastSuccess", "PastDefection", "Region(Americas)", "Region(East Asia + Pacific)", "Region(Europe + Eurasia)", 
                             "Region(Middle East + North Africa)", "Region(South + Central Asia)")

stargazer(DurationTable, title = "Surival Analysis",
          align=TRUE, float = TRUE, column.sep.width = "0pt")
